\chapter{Segmentaci\'on de im\'agenes}

\section{{\color{blue}Introducci\'on}}
\label{sec:introduccion}

La segmentaci\'on de las diferentes regiones del cerebro a partir de im\'agenes
de resonancia magn\'etica es un proceso clave a la hora de realizar estudios
sobre pacientes, tanto en casos sanos como aquellos con la presencia de tumores.
Automatizar este proceso es una tarea complicada y su validaci\'on suele
realizarse mediante la comparaci\'on de los resultados arrojados por el algoritmo
y los obtenidos manualmente; estos \'ultimos son realizados por expertos y la
forma de obtenci\'on resulta costosa dada a una combinaci\'on de factores que se
encuentran presentes en la mayor\'ia de los casos. Algunos de estos factores
resultan ser:

\begin{itemize} 
 \item La deformaci\'on de tejidos sanos debido a la presencia de masa
tumoral reduce dr\'asticamente la posibilidad de utilizar la informaci\'on
obtenida de estudios sin la presencia de tumor o edema.
 \item Las transiciones graduales en las intensidades de un tumor o edema a un
tejido sano afectan la precisi\'on de las segmentaciones, ya sea utilizando modelos de detecci\'on de
bordes, de regiones o realizando el proceso manualmente.
 \item Las diferentes modalidades normalmente
 utilizadas en estudios de resonancia magn\'etica no son capaces de
 diferenciar por si solas todos los tejidos (las im\'agenes T1 reforzadas con
 gadolinio suelen devolver intensidades de se\~nal similares tanto en la presencia de l\'iquido 
 cefalorraqu\'ideo como en la de sangre, y ciertas partes del tumor que
 pertenecen al tejido necr\'otico no aparecen resaltadas en esta modalidad).
\end{itemize}

%\cite{Cobzas2007}

\section{{\color{blue}Estado del arte}}
\label{sec:stateOfArt}
% comentar un poco los metodos y que no tienen y que cosa tiene el nuestro que
% estos no.

% seguir introduccion de a partir de \cite{Boykov2006}

En los \'ultimos veinte a\~nos, en el \'area de visi\'on en computaci\'on, se han
implementado con \'exito una gran cantidad de nuevos algoritmos con el objetivo
de segmentar diferentes objetos en im\'agenes. Modelos como el de snake
\cite{Kass1988}, \cite{Cohen1991}, active contours \cite{Isard1998}, geodesic
active contour \cite{Caselles1997}, \cite{Yezzi1997} y t\'ecnicas de ``shortest
path'' \cite{Mortensen1998}, \cite{Falcao1998} son solo un peque\~no grupo
representativo de tantos otros que buscan particionar la imagen de acuerdo a
determinados objetivos, partiendo de un conjunto de caracter\'isticas. La
t\'ecnica que presentamos en este estudio se basa en la minimizaci\'on de
energ\'ia mediante graph cuts y max-flow que ya ha demostrado otorgar resultados
prometedores en problemas de visi\'on general. Ejemplos de esto es la
restauraci\'on de im\'agenes \cite{Greig1989}, \cite{Boykov2001a}, stereo
\cite{Roy1998}, \cite{Ishikawa1998}, \cite{Boykov1998}, multi-view reconstruction
\cite{Kolmogorov2002}, \cite{Vogiatzis2005}, \cite{Lempitsky2006} o texture
s\'intesis \cite{Kwatra2003}. Nuestro estudio apunta a incorporar la
informaci\'on contenida en im\'agenes multi-modales al algoritmo de
segmentaci\'on mediante graph cuts y la minimizaci\'on de energ\'ia a partir del
modelo de Mumford-Shah. La efectividad de este tipo de algoritmos ya fue
demostrada en una gran cantidad de estudios recientes sobre problemas de visi\'on
basados en los conceptos que plantean Boykov y Jolly \cite{Boykov2001a}. Los
m\'etodos de segmentaci\'on basados en el concepto de energ\'ia y su
minimizaci\'on pueden distinguirse por el tipo de funci\'on que utilizan para
realizar esta \'ultima. Una divisi\'on general plantea dos grupos; por un lado
los m\'etodos que utilizan un funcional basado en un contorno continuo y por
otro, aquellos que trabajan con datos discretos. Nuestro modelo pertenece al
segundo grupo y se caracteriza por manejar un conjunto finito de variables
asoci\'andolas con los nodos del grafo que representan un voxel de la imagen a
segmentar. Al trabajar con este tipo de modelos discretos y gracias a la
utilizaci\'on de graph cuts, somos capaces de incorporar im\'agenes en 3D sin
necesidad de modificar la forma en que trabaja el algoritmo.\\ Previo a la
utilizaci\'on de graph cuts para la segmentaci\'on de im\'agenes, la
optimizaci\'on mediante la busqueda de m\'aximos$-$m\'inimos globales era
posible solo para algunos m\'etodos de segmentaci\'on en 2D. En general, las
soluciones globales resultan ser m\'as atractivas debido a su mayor estabilidad,
potencialmente hablando. Por ejemplo, las imperfecciones en soluciones obtenidas
mediante \'optimos globales est\'an relacionadas al coste de la funci\'on en
lugar de ser derivadas de problemas num\'ericos durante la minimizaci\'on
\cite{Boykov2006}. Algunos algoritmos como los de active contours
\cite{Cohen1997}, shortest path \cite{Mortensen1998}, \cite{Falcao1998}, ratio
regiones \cite{Cox1996}, y algunos otros m\'etodos \cite{Jermyn1999} computan una
soluci\'on mediante \'optimos globales en una aplicaci\'on 2D cuando los bordes
de la segmentaci\'on son una curva de 1D. Y dado que ninguno de estos m\'etodos
son capaces de ser aplicados sobre segmentaciones en 3D, las posibles soluciones
implican la utilizaci\'on de m\'etodos ad-hoc (segmentando los cortes uno por
uno) o de t\'ecnicas que aproximen mediante \'optimos locales. Nuestra
implementaci\'on se encuentra dentro del segundo grupo bas\'andonos en las
propiedades que describe \cite{Boykov2001} sobre los m\'inimos locales obtenidos
de determinada forma.\\
Num\'ericamente hablando, nuestra implementaci\'on se basa en un algoritmo de
graph cuts sin modificaciones utilizando una optimizaci\'on combinatoria sobre un
modelo de energ\'ia discreta \cite{Ford1962}, \cite{Goldberg1988},
\cite{Boykov2004}. Y a diferencia de otros m\'etodos que se basan en esquemas de
aproximaci\'on num\'erica y que deben ser cuidadosos en la etapa de dise\~no para
poder asegurar un alto nivel de calidad, la utilizaci\'on del algoritmo de graph
cuts sin modificaci\'on nos asegura el necesario nivel de calidad sin mayores
inconvenientes.\\
El estudio de \cite{Boykov2004} muestra la eficiencia de los algoritmos de graph
cuts y max-flow aplicados a problemas de visi\'on. A su vez se ha demostrado que
las t\'ecnicas de max-flow pueden resolver problemas de segmentaci\'on en 2D y 3D
en aproximadamente tiempo real usando computadoras de escritorio. Aunque este
tipo de implementaciones requieren una gran cantidad de memoria, resultados presentados
por \cite{Lombaert2005} muestran que determinadas t\'ecnicas (multilevel y
banded) pueden aliviar el problema.\\
Recientemente \cite{Bagon2006} mostr\'o una relaci\'on directa entre graph-cuts y
level-sets \cite{Sethian1999}, \cite{Osher2002}, \cite{Sapiro2001},
\cite{Osher2003}; los resultados presentados en el estudio muestran que los
algoritmos basados en graph cuts pueden ser usados en los problemas frente a los
que antes se aplicaban m\'etodos basados en level-set. Dado que la utilizaci\'on
de graph cuts permite optimizar formulaciones discretas de energ\'ia combinando
regularizaci\'on de bordes y de regiones de manera similar a la planteada por el
funcional de Mumford-Sha en el area de problemas continuos, los problemas de
integraci\'on de regiones, detecci\'on de bordes y formas pueden ser resueltos
mediante este tipo de implementaciones sin mayores inconvenientes. Dado que
esta implementaci\'on se basa en una representaci\'on implicita
de los bordes de los objetos en la imagen, las segmentaciones obtenidas poseen
las mismas propiedades topol\'ogicas que aquellas resultantes de la utilizaci\'on
de level-set.

\section{{\color{blue}Funcional de Mumford-Shah}}

Una de las principales dificultades en el area de visi\'on en computaci\'on es
entender como se realiza la conversi\'on de la se\~nal a un observable
geom\'etrico. Sea $I(x)$ (donde x representa la tupla $(x,y)$) una imagen
definida en el dominio $\W$ y sin ning\'un tipo de estructura geom\'etrica a
priori; nuestro objetivo consiste en particionar a $\W$ en diferentes $\W_i$.
Siendo estas \'ultimas delimitadas por un borde regular que llamaremos
$\K$, y en las cuales la se\~nal de $I$ debe ser homogenea.\\
Se han desarrollado una gran cantidad de algoritmos de segmentaci\'on que buscan
particionar la imagen en regiones homogeneas delimitadas por bordes continuos.
Esta tarea no resulta trivial debido a que las regiones y los bordes son
entidades de dimensiones diferentes, geom\'etricamente hablando (las regiones
poseen dos dimensiones mientras que los bordes solo una). Pero detr\'as de esta
cantidad de estudios se suele encontrar variaciones del mismo modelo. Seg\'un
\cite{Morel1995} la mayor\'ia de los algoritmos de
segmentaci\'on trabajan minimizando la misma energ\'ia.\\
La idea de energ\'ia nace de la necesidad de comparar una segmentaci\'on con
otra, y de medir cuan bien aproximan a la imagen original $I$. Uno de los modelos
de energ\'ia m\'as conocidos, y sobre el cual nos basamos, es el de Mumford-Shah.
En \cite{Mumford1989} se plantea que en el modelo Bayesiano existen dos partes:
el modelo a priori y el de datos. El modelo a priori se basa en el conocimiento
de lo que se define como segmentaci\'on, que resulta ser una aproximaci\'on de la
imagen original $I$ mediante funciones suaves a pedazos $u$ en $\W - \K$ y
discontinuas en el conjunto de bordes $\K$. Nuestra tarea consiste en elegir de
todas las aproximaciones $(u,\K)$ de $I$, la mejor posible. Para esto
\cite{Mumford1989} define un funcional de energ\'ia $E(u,\K)$ que posee tres
t\'erminos:

\begin{enumerate}
\item Un primer t\'ermino que controla que $u$ sea suave en las regiones $\W_i$
de $\W - \K$.
\item Un segundo t\'ermino que controla cuan bien aproxima $u$ a la imagen.
\item Y un \'ultimo t\'ermino que controla la longitud y la suavidad de $\K$. 
\end{enumerate}

El funcional entonces resulta:

\begin{equation} \label{eq:msEnergy}
E(u,\K) = \int_{\W-\K} | \nabla u |^2 dx + \lambda \int_{\W} (u - I)^2 + \mu \int_{\K} d\sigma
\end{equation}

Mediante la variaci\'on de los coeficientes $\lambda$ y $\mu$ logramos diferentes
tipos de segmentaciones. Por ejemplo a mayor $\mu$ disminuiremos
la granularidad de la segmentaci\'on.\\
Como podemos ver, la minimizaci\'on de $E$ implica lo siguiente:
\begin{enumerate}
\item A medida que se minimiza la energ'ia, $u$ resulta ser mas homogenea en las
regiones $\W_r$. Si por ejemplo $u$ es constante entonces $\nabla u = 0$ y 
$\int_{\W-\K} | \nabla u |^2 dx = 0$, de modo que minimizando el primer t\'ermino
 obligamos a $u$ a ser lo mas constante posible en los $\W_r$.
\item Al mismo tiempo se obliga a que $u$ aproxime lo suficientemente ``bien'' a
$I$. Por ejemplo si $I = u$ entonces $\int_{\W} (u - I)^2 = 0$, de modo que minimizar el
segundo t\'ermino fuerza a $u$ a parecerse lo m\'as posible a $I$.
\item Finalmente, mediante $\K$ se mide la longitud de los bordes, y a medida
que se minimiza la energ\'ia se se reduce la probabilidad de
sobre-segmentaci\'on.
\end{enumerate}

A partir de este modelo gen\'erico se han desarrollado diferentes
aproximaciones utilizando funciones cuyo comportamiento responde al descripto en
\cite{Mumford1989}. Uno de ellas y sobre la cual nos basaremos en este estudio,
es la presentada en \cite{Petitot2003}:

\subsection{{\color{blue}Formulaci\'on probabil\'istica}}
\label{sec:probability}

La formulaci\'on mediante probabilidades se basa en asumir que las intensidades
de los v\'oxeles dentro de cada region resultan de la aplicaci\'on de un proceso
aleatorio dada una funci\'on de densidad. Sea $p_i$ la distribuci\'on de
probabilidad en la regi\'on $\Omega_i$, en el estudio realizado por
\cite{Paragios2002} se muestra que al asumir que
las intensidades de los v\'oxeles responden a un proceso aleatorio, la maximizaci\'on de la probabilidad de una
partici\'on de I se logra mediante la minimizaci\'on de la energ\'ia:

\begin{equation}
E(\{\Omega_i,i= 1..N\}) = - \sum_i \int_{\Omega_i} \log p_i(I(x)) dx
\end{equation}
		
La estructura de la imagen solo se logra observar mediante las regiones $i$, y el
funcional no posee informaci\'on para limitar las caracter\'isticas espaciales.
Agreg\'andole un t\'ermino regularizacion para el borde C entre las regiones
resulta:

\begin{equation} \label{eq:probabilidad}
E(\{\Omega_i,i= 1..N\}) = - \sum_i \int_{\Omega_i} \log p_i(I(x)) dx + \nu|C|
\end{equation}

Para llegar a expresar el funcional de esta manera necesitamos hacer una serie de
asumpciones. A partir del estudio realizado por \cite{Paragios2002}, consideramos a
$p(P(\Omega)|I)$ como la probabilidad a posteriori dada la imagen $I$. Luego, la
partici\'on \'optima de la imagen se obtiene mediante la maximizaci\'on de esta
probabilidad, y mediante la regla de Bayes podemos expresar a esta de la
siguiente forma:

\begin{equation} \label{eq:bayes_1}
p(P(\Omega)|I) = \frac{p(I|P(\Omega))}{p(I)}p(P(\Omega))
\end{equation}

siendo $p(P(\Omega))$ la probabilidad de la partici\'on $P(\Omega)$ y $p(I)$ la
probabilidad de la imagen $I$. El tercer t\'ermino, $p(I|P(\Omega))$ representa
la probabilidad a posteriori de la segmentaci\'on de la imagen $I$, dada la
partici\'on $P(\Omega)$. Al segmentar una imagen, el t\'ermino  $p(I)$ se
vuelve constante y \mathref{eq:bayes_1} resulta:

\begin{equation} \label{eq:bayes_2}
p(P(\Omega)|I) \propto p(I|P(\Omega)) p(P(\Omega))
\end{equation}

La expresi\'on \ref{eq:probabilidad} es la base de varios estudios
\cite{Leclerc1989}, \cite{Zhu1996}, \cite{Samson2000} y \cite{Paragios2002}.
Luego de las aproximaciones realizadas para obtener esta expresi\'on, es importante
tener presente que condiciones debe satisfacer la imagen para poder ser
segmentada mediante este enfoque. Las principales simplificaciones y sus
consecuencias:

\begin{itemize}
  \item Los bordes de las particiones se manejan mediante una regularizaci\'on
  de los mismos.
  \item El dominio $\Omega$ de la imagen se compone de $N$ regiones, y este
  n\'umero debe ser conocido \textit{a priori}. Esto l\'imita fuertemente la
  posibilidad de crear un proceso desatendido.
  \item Las intensidades de los p\'ixeles dentro de una misma regi\'on son el
  resultados de un proceso aleatorio, identicamente distribuidos e
  independientes entre si. Esta es la asumci\'on m\'as fuerte de todas y
  restringe el conjunto de imagenes a aquellas que se aproximan a ser
  constantes de a pedazos.
\end{itemize}

\subsubsection{{\color{blue}Los par\'ametros de la funci\'on de densidad}}

La definici\'on de un criterio gen\'erico sin plantear ning\'un tipo de
hip\'otesis sobre la distribuci\'on de las intensidades de los p\'ixeles en una
regi\'on es imposible. Frente a esto existen varias posibles decisiones a
tomar. Una forma sencilla de atacar esta necesidad de caracterizar las
distribuciones es la de asumir que todas responden a una misma funci\'on de
densidad pero con diferentes par\'ametros. Luego el criterio de segmentaci\'on
se puede optimizar con respecto a dos conjuntos de parametros: la partici\'on
de la imagen y los par\'ametros estad\'isticos. A partir de la f\'ormula
\mathref{eq:bayes_3} podemos definir el criterio de optimizaci\'on
introduciendo los par\'ametros estad\'isticos como desconocidos. Si
$p(I|\theta_i)$ es la representaci\'on de la probabilidad a posterior en la
regi\'on $\Omega_i$ parametrizada por $\theta_i$, entonces la segmentaci\'on
buscada se obtiene de la minimizaci\'on de:

\begin{equation}
E(\{\Omega_i, \ldots,\Omega_N\},\{\theta_i, \ldots,\theta_N\}) = - \sum_i
\int_{\Omega_i} \log p_i(I(x)|\theta_i) dx + \nu|C|
\end{equation}

Esta energ\'ia tiene dos tipos de parametros: las particiones de la imagen
$\Omega_i$ que son parte de $\Omega$, y los par\'ametros estad\'istico
$\theta_i$ que pertenecen al espacio de par\'ametros $\Theta$.

\subsection{{\color{blue}Aproximaci\'on mediante Gaussianas multivaluadas}}

Una forma directa al trabajar con varias modalidades de una im\'agen es la de
utilizar como aproximaci\'on a una funci\'on de densidad Gaussiana
multivaluada:

\begin{equation}
p(z|\mu, \Sigma) =
\frac{1}{(2\pi)^{d/2}|\Sigma|}e^{-\frac{1}{2}(z-\mu)^T\Sigma^{-1}(z-\mu)}
\end{equation}

Como dijimos en \ref{sec:probability}, la aproximaci\'on mediante esta funci\'on
de densidad limita a que las intensidades de los p\'ixeles dentro de una regi\'on
se distribuyan de manera normal. En \cite{Fischl2002} se presenta un estudio
sobre la distribuciones de intensidades de los v\'oxeles en los diferentes
tejidos del cerebro. El mismo muestra que la mayor\'ia de las regiones
se comportan de manera aproximada a la que se obteniene mediante una funci\'on
de densidad Gaussiana \viewref{fig:distribucion} \footnote{\cite{Fischl2002}}.

\insertImage{im/distribucion}{Histograma de intensidades para materia blanca
(WM), materia gris (GM), ventriculo lateral (lV), t\'alamo (Th), Caudate (Ca),
Putamen (Pu), Pallidum (Pa), Hippocampus (Hp), and Amygdala
(Am)}{fig:distribucion}{100mm}
 

